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The statistical significance of network properties is conditioned on null models which satisfy spec¬ 
ified properties but that are otherwise random. Exponential random graph models are a principled 
theoretical framework to generate such constrained ensembles, but which often fail in practice, ei¬ 
ther due to model inconsistency, or due to the impossibility to sample networks from them. These 
problems affect the important case of networks with prescribed clustering coefficient or number of 
small connected subgraphs (motifs). In this Letter we use the Wang-Landau method to obtain a 
multicanonical sampling that overcomes both these problems. We sample, in polynomial time, net¬ 
works with arbitrary degree sequences from ensembles with imposed motifs counts. Applying this 
method to social networks, we investigate the relation between transitivity and homophily, and we 
quantify the correlation between different types of motifs, finding that single motifs can explain up 
to 60 % of the variation of motif profiles. 

PACS numbers: 05.10.Ln, 64.60.aq, 89.75.He 


Networks form the basis of an ample class of com¬ 
plex systems. The observed topological patterns of such 
systems often yield the only available evidence for the 
underlying principles behind their formation. However, 
the significance of any observed property can only be 
assessed in comparison to a properly defined network en¬ 
semble that acts as a “null” model m ■ For instance, 
clustering (i.e. high density of triangles), skewed degree 
distributions, and community structure are considered 
significant in real networks because they are absent in 
Erdos-Renyi networks. To perform such comparisons, it 
is essential not only to properly define such null mod¬ 
els, but also to correctly sample network realizations 
from them. This is relatively straightforward when the 
ensemble generates networks where the edges are sam¬ 
pled independently (e.g. Erdos-Renyi and configuration 
models HE], the stochastic block model mm) and it 
remains feasible when strict edge independence is vio¬ 
lated due to hard constraints [844TO]. However, for en¬ 
sembles with more generic constraints the sampling is 
significantly more challenging. A particularly important 
example is ensembles with a prescribed density of con¬ 
nected subgraphs (“motifs”) mm- For this class of 
models, one often finds abrupt phase transitions, where 
sampled networks possess either very high or very low 
motif density mm , excluding intermediary values of¬ 
ten encountered in real systems. Furthermore, they often 
show strong non-ergodic behavior, with very slow relax¬ 
ation that forbids unbiased sampling in practical compu¬ 
tational time m- Since the edge placement is not inde¬ 
pendent, the densities of different motifs are correlated 
with each other and also with large-scale network struc¬ 
tures nn nsi- Without addressing the issue of correct 
sampling, these correlations cannot be properly identi¬ 
fied, which makes the occurrence of these patterns in 
real systems difficult to interpret. In particular, it is 
not possible to conclude whether a particular motif den¬ 


sity profile indicates a topology optimized towards ro¬ 
bustness Hi 17] or whether it is merely a byproduct of 
a specific large-scale structure [Hi nsi, of combinatorial 
constraints [TO], or of correlations between motifs. 

In this Letter we show how to sample from ensembles 
with prescribed motif densities in polynomial time. We 
employ a multicanonical Monte Carlo method |2J that 
allows the entire range of the order parameter to be 
explored. In this manner, not only the non-ergodicity 
problem is explicitly avoided, but it also becomes pos¬ 
sible to sample networks with arbitrary motif densities, 
even those at intermediate values that are unattainable 
via traditional importance sampling. This allows us to 
quantitatively investigate two fundamental problems in 
social networks: the homophily-transitivity relationship 
and the interdependence of different motif types. 

We are interested in network ensembles that possess 
one particular observable s of interest, but that are oth¬ 
erwise maximally random. The last requirement is es¬ 
sential to ensure that the ensemble is representative of 
the networks with a given s and is not subject to addi¬ 
tional (hidden) constraints. Both features are achieved 
by sampling the network from an exponential random 
graph model (ERGM) [31 [10l [221425] Q , where each graph 
g £ Q occurs with probability 

„/3s(g) 

n /3 (g) = ~ , where Z 0 = V e 0s( - 9 \ (1) 

Zp 9 es 

where s(g) is the observable associated with network g , 
and /? is an inverse-temperature parameter, in analogy 
to the canonical ensemble in statistical physics. The 
distribution of s is pp(s) = J2 g egfi( s (g) ~ s)Up(g) = 
po(s)eP s /Z 0 , where po(s) = pp=o(s) is called the state 
density (the fraction of networks g in the ensemble that 
have observable equal to s). The ensemble that acts as a 
null model for an empirical network with s = s* is usually 
constructed fixing /3 in such a way that (s)p = ]C S sp 0 (s) 
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FIG. 1. (Color online) Multicanonical sampling of exponen¬ 
tial random graphs with imposed clustering avoids the limi¬ 
tations of canonical sampling. The ensemble Q corresponds 
to k-regular undirected networks with N = 640 and degree 
k = 4. The observable is the clustering coefficient s = c 
(proportional to the number of triangles, n a). The main 
plot shows (c) (and standard deviation) as a function of the 
inverse temperature fl obtained by canonical (symbols) and 
multicanonical (continuous thick line) sampling. Inset: the 
distribution pp{c ) for /3 = 3.54 ss /3 pt obtained by the two 
methods. Canonical samplings used 5 x 10 5 MCMC steps for 
equilibration, before another 5 x 10 s steps were used for esti¬ 
mation. After these steps, the value of /? was slowly increased 
(/3 t) or decreased (/3 4-) and the process repeated. The mul¬ 
ticanonical sampling used 20 Wang-Landau steps to estimate 
Po(c) (each step used 5 tunneling steps) [21. 121). 


equals s*. The number of networks in this ensemble 
typically grows exponentially with the number of nodes, 
and, thus, besides a small set of observables s that can 
be treated analytically, investigation of ERGMs requires 
sampling networks g from Q using Monte Carlo meth¬ 
ods 123], 

The usual approach of sampling from Q is via Markov 
chain Monte Carlo (MCMC) method works as follows: 
starting from one network g £ G, a new network g' £ Q n 
is proposed by choosing two links at random and ex¬ 
changing one of the nodes of each link, which preserves 
the degree-sequence of the network [26]. The proposed 
network is accepted with the Metropolis-Hastings prob¬ 
ability A(g i —y g') = min{l, e^ a ^ 9 ) -s (sb} and the pro¬ 
cess is repeated from g' ( g ) if the proposal is accepted 
(rejected) [22] - Since the moves fulfill ergodicity and de¬ 
tailed balance, for sufficiently long times the values of 
s in the sampled networks g are distributed as pp(s). 
However, despite this asymptotic guarantee, in practice 
this method often fails because the time to approximate 
pp (s) grows exponentially with the number of nodes N. 
This happens whenever pp possesses more than one local 
maximum (minimum of the free energy) and the barri¬ 


ers between them grow with N. As we show below, this 
generically happens when the observables s are related 
to motifs. 

As an alternative to the canonical (simple Metropo¬ 
lis) sampling method described above, we propose a 
multicanonical sampling to overcome the aforementioned 
problem. This method aims to sample networks uni¬ 
formly on a pre-defined observable range [s ra m,s mM ], 
thus overcoming the minima of pp(s) that are respon¬ 
sible for the weak performance of the canonical method. 
This is done by sampling the states according to aux¬ 
iliary ensemble with probabilities n'(g) oc l/po(s(g)), 
achieved by simply changing the acceptance to A(g H > 
g') = min{l, po(s(g'))/ Po(s(g))} 0- However, in order 
to perform this sampling we need to know the state den¬ 
sity po{s)- In order to estimate it, we use the Wang- 
Landau algorithm mm, which, in short, constructs an 
adaptive histogram to approximate po(s(g)) 0- After 
convergence, pp{s) is estimated for all /3’s reweighting 
Po(s) through pp (s) = po(s) exp(— (3s)/Zp 0. Hence, the 
auxiliary ensemble allows to explore the original canoni¬ 
cal ensembles without being restricted to the most prob¬ 
able regions. More importantly, we can impose the de¬ 
sired value of the observable as a hard constraint a pos¬ 
teriori, i.e., only sample networks with s(g) = s*. The 
multicanonical approach has recently been applied to in¬ 
vestigate the spectral gap of networks 03], and related 
approaches have been used to investigate percolation [31] 
and resilience properties of networks [32 . 

In Fig. [l] we show how the application of multicanoni¬ 
cal sampling solves the limitations of canonical sampling 
in the classical problem of introducing clustering in a k- 
regular network mm- Here, nodes are forced to have 
the same degree k and the observable of interest is the 
number of triangles, s(g) = ha- Fixing ua is the same 
as fixing the clustering coefficient c = 3ua /^a > where tia 
is the number of connected triples (a constant for all net¬ 
works with the same degree sequence) [3] . This model ex¬ 
hibits a transition at a specific value of ft = f3pr (~ 3.54 
for k = 4), separating low and high-clustering phases [T3] . 
The canonical sampling is unable to compute (c) close 
to the phase-transition because it yields different esti¬ 
mations of (c), depending whether f3 is slowly increased 
(/3 ti lower branch) or decreased (/? 4-! upper branch). 
This hysteresis is typical around first-order phase transi¬ 
tions (coexisting phases) and indicates that the canonical 
sampling is in a metastable state. Indeed, pp PT {c) has 
two local maxima in which the canonical sampling be¬ 
comes trapped (inset in Fig. [l]). On the other hand, the 
multicanonical sampling is immune to these problems: it 
correctly characterizes (c) at /3 = /?pt and reveals the full 
distribution pp-p PT . Hence, the method is not only ca¬ 
pable of computing the correct ensemble average for any 
(3, it yields typical networks with any value of c, including 
the significant gap c £ [.2,1] which is unattainable with 
the canonical sampling. In Fig. [2] we confirm that the 
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FIG. 2. (Color online) Efficiency of the multicanonical 
method to sample networks with constraints. The compu¬ 
tational cost (in number of MCMC steps) to generate an in¬ 
dependent realization of a network in the fc-regular ensemble 
with k = 3 is plotted as a function of N. In the canoni¬ 
cal method close to the critical /?, this requires passing the 
minimum of pp(c) (inset of Fig. [l]). We measured that the 
height of this barrier increases as A p sa 0.4./V, which leads 
to an exponential increase in the cost (dashed line). Sam¬ 
pling independent realizations in the multicanonical method 
requires, at most, a tunneling (the number of MCMC steps 
todoc = 0^>c=l^>c = 0) [213 - The measured tunneling 
time (circles and full line) scales polynomially. Inset: con¬ 
vergence of the relative error in the logarithm of the density 
of states (entropy) during convergence of the Wang-Landau 
algorithm, estimated comparing the measured value with the 
exact value on c=l. The saturation of the error observed for 
large number of steps does not hinder the sampling of any c 
(see Ref. [33] for methods to overcome the saturation). 

computational cost of the multicanonical method scales 
polynomially with system size, a dramatic improvement 
over the exponential scaling of the canonical method. 

Next we use the multicanonical method to investi¬ 
gate two important problems of social networks. The 
first problem we consider is to distinguish between ho- 
mophily (the tendency of “similar” nodes to connect to 
each other) and transitivity (the tendency of nodes that 
already share a common neighbor to connect to each 
other) in social networks 2] HU 13HH38] , We use the 
(undirected) network of email exchange within a univer¬ 
sity [M]. It consists of N = 1,133 users, and M = 5,451 
email exchanges, and a roughly exponential degree dis¬ 
tribution. As observables we consider the clustering co¬ 
efficient c and the degree assortativity r [3fJ], for which 
we obtain c* = 0.166(12) and r* = 0.08(3) (uncertainties 
in the last digit estimated using the order-10 Jackknife 
method). We assess the significance of these values by 
comparing them to those obtained in the following three 
network ensembles with the same degree sequence as in 
the original network: 

(i) Same weight to all networks g (i.e. the configuration 
model). Canonical sampling with j3 = 0 yields (c)^—q = 
0.028(1) and {rp— o) = —0.017(13), much smaller than c* 
and r* as typically found in social networks. 



FIG. 3. (Color online) Relationship between clustering c 
and assortativity r in the email network of Ref. [31]. The 
assortativity r* is higher than in a random network with the 
same degree sequence, but lower than in a typical network 
with a fixed clustering c*. The plot shows r and c of different 
networks: the email network (green), a typical fully random 
network (red), a typical random network with c = c* (black), 
and networks sampled using the multicanonical method from 
an ensemble with equal probability for networks with the same 
c (blue dots). Inset: (c) obtained using the canonical method. 


(ii) ERGMs with (c) = c*. In order to determine whether 
the assortativity is a consequence of high clustering tnj 
we would like to measure (r) from the null model with 
(c) = c*. This canonical sampling fails because (cb vs. 
(3 shows an hysteresis around s = c* (inset of Fig. [2j in 
agreement with our previous discussion). 

(iii) Hard constraints with c(g) = c *, obtained using mul¬ 
ticanonical sampling. As mentioned before, this type 
of hard constraint is unfeasible with canonical sampling, 
even if the desired observable value is realizable. With 
the multicanonical method we sample points after a num¬ 
ber of Monte Carlo steps proportional to the tunneling 
time, which guarantees that the sampled points are inde¬ 
pendent and unbiased [21j . We performed multicanonical 
sampling for a desired c and measured the assortativity 
r. The results are shown in Fig. [3] and reveal that ran¬ 
dom networks with the same clustering of the email net¬ 
work c = c* typically show a much larger assortativity 
(r) > r*. Therefore, although both c* and r* are larger 
than one would expect for a fully random network, the 
actual value of r* is significantly less than one would ex¬ 
pect by knowing only c*. From this we conclude that the 
degree homophily is not explained alone by transitivity. 

The second problem we address is the extent to which 
the occurrence of different motifs (connected subgraphs) 
are related to each other and the impact of such correla¬ 
tions on the so-called motif profiles El- Here we focus on 
directed networks, and the observable of interest is the 
number rii of occurrences of a specific motif i. Again, 
traditional sampling methods are not suited to address 
this problem because of the existence of (potentially mul- 
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FIG. 4. (Color online) Motifs are correlated to each other in blocks, (a) The Pearson correlation coefficient [Rij = (( rurij) — 
{ni){rij))/an i crn j ] between motifs i and j , computed by varying the constrained motif (in a range which includes the values of the 
real and random networks seem [TT] ). (b) Upper panel: Motif profile [T7] built from the z-score Zj vs. j (Zj = ( rij — {nj))/a nj , 
where rij is the number of motifs j and (...) and ay are the average and standard deviation in the /3 = 0 ensemble). Different 
lines correspond to the Zj of the real network (z], blue line) and the expected Zj’s in the constrained ensemble in which m is 
equal to the n* of the real network, where i is the constrained motif shown in the legend (zj = Zj for j = i). Middle panel: the 
correlation between the profiles shown in the upper panel, i.e., between the profile 2 of the real network and the profile z' of the 
motif-i constrained network (as a function of i), computed as R zz / = (( zz') — (z)(z'))/a z , where (...) and o z were computed 
over j ^ i. Lower panel: comparison of the z-score shown in the upper panel (blue line) and the alternative z-score obtained 
computing (...) and oy in the ensemble constrained by rii = n*, where i is indicated in the legend (Zj = 0 for j = i). 


tiple p)3] ) discontinuous phase transitions. Instead, using 
the multicanonical method, we reliably sample networks 
with a prescribed count of one particular motif. By mea¬ 
suring the counts of all other motifs, we obtain the corre¬ 
lations between them and the constrained motif. In this 
manner, we obtain |41j the interdependence between all 
13 different 3-node motifs in a directed acquaintance net¬ 
work between physicians 3D] (with N = 241 nodes and 
M = 1,098 edges). The results in Fig. [4] reveal strong 
positive and negative correlations between pairs of mo¬ 
tifs. Two blocks of motifs can be identified (1,2, 3, 7, 8 
and 4,5,6, Fig. |4ji). Motifs show positive correlations 
within their blocks and are anti-correlated with motifs 
in the other blocks (the motifs 9 to 13 show a mixed 
behavior). Given that one motif is over (under) repre¬ 
sented, one should expect also an over representation in 
motifs positively (negatively) correlated with it. As a 
consequence of this correlation, we find that single mo¬ 
tifs explain up to 60% of the variance of the motif profile 
across the other 12 motifs (Fig. upper and middle 
panel). Furthermore, if the constrained ensembles are 
used to compute alternative z-scores, we find that the re¬ 
sulting motif profiles vary dramatically depending on the 
constraint, with some motifs j showing variations from 
Zj » 0 to Zj « 0 (Fig. lower panel). This sensitiv¬ 
ity of the motif profile Zj shows that such profiles bring 
limited insights on the over- or under-representation of 
individual motifs in a network. In particular, since such 
non-trivial profiles as those seen in Fig. & can be ob¬ 


tained by imposing the occurrence of a single motif, it is 
questionable whether conclusions regarding the underly¬ 
ing formation mechanisms can be reliably reached from 
them HU EH]. Nevertheless, the null models considered 
here represent a principled approach of assessing the rela¬ 
tive significance of motif occurrences that is more mean¬ 
ingful than the usual comparison to fully random net¬ 
works. 

In summary, we have shown that multicanonical sam¬ 
pling allows for an improved network generation and for 
the investigation of problems which were otherwise in¬ 
tractable. In particular, we characterize ERGMs in cases 
where the usual canonical sampling fails and we sam¬ 
ple networks imposing hard constraints, an alternative 
to a direct sampling of ERGMs even when the usual 
algorithms are feasible. Our analysis of empirical net¬ 
works demonstrates that using the multicanonical sam¬ 
pling allows the investigation of the interdependence be¬ 
tween network properties. In particular, we quantified 
the correlation between clustering and assortativity, and 
between different motifs, as well as the extent to which 
their significance profiles can be explained by single mo¬ 
tifs. This opens the possibility of investigating the corre¬ 
lation between motifs as well as other local-scale proper¬ 
ties and the large-scale structure of networks [14] , such as 
communities, core-peripheries and many others. The sys¬ 
tematic disentangling of these diverse features is a crucial 
and open problem in the identification of fundamental 
models of network formation. 
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Supplemental Material 

Wang-Landau algorithm to sample networks 

The sampling algorithms used in the paper perform 
random walks in the space of constrained networks. In 
our case, this space is built by networks with a fixed num¬ 
ber of nodes and a fixed degree sequence (e.g. the degree 
sequence of the e-mail network). Besides this space, the 
observable we are interested in characterizing s(g) (e.g. 
number of triangles of the network g ) and its range of 
interest [s m i n , s max ] are also chosen a priori. 

The Wang-Landau algorithm performs a random walk 
in the space of constrained networks that aims to visit 
equally often any value of s £ [s m i„, s max ]. The outputs 
of this algorithm are: 1. a numerical approximation of 
the state density po(s); and 2. a set of random networks g 
such that their observables s(g) are uniformly distributed 
in [s m in,s max ]. Below we describe the main steps of the 
Wang-Landau algorithm. 

The following quantities have to be initialized and 
evolve in time: 

• a network g, initially set to be an arbitrary network 
with s(g) e [s min ,s max ]; s = s(g) represents its 
observable. 

• a histogram-like list S(s), for s £ [s m in, s m ax] 
(binned if s is continuous), that represents a dis¬ 
crete approximation of logp(s); it is initialized for 
all s to S(s) = 0. 

• the Wang-Landau refinement parameter /, initial¬ 
ized at / = 1 (the minimum value / m i n is set a 
priori, e.g. / min = 2“ 12 ). 

The algorithm evolves according to the following rules: 

1. Repeat until a pre-defined number (e.g., 10) of 
round-trips are achieved: 

(a) Randomly propose a new network g' in the 
space of constrained networks, and compute 
s' = s(g’) (see how below); 

(b) Update g, s to g = g', s = s' if log(r) < 
a(s',s) = S(s') — S(s) where r is a random 
number drawn from a uniform distributed in 

[o,i]; 

(c) Update S(s) = S(s) + f. 

2. update / to / = //2, go to 1. if / > / mi „. 

After the convergence of the evolution (/ < / m i n ) de¬ 
scribed above, S(s) is an approximation of logp(s) up to 
a normalization constant. Setting / = 0 at this point and 
repeating loop 1 generates networks g such that s(g) is 
uniformly distributed in [s m i n ,s max ] (multicanonical en¬ 
semble) . 

Notes: 


i) A round-trip is achieved in step 1. when, for the 
first time, the observable went from s = s m j n to 
s = Smax and returned back. 

ii) The network g' is constructed by selecting two 
edges of the original network g and exchanging one 
of the nodes of one edge by one of the nodes of 
the other edge. This guarantees that the degree 
sequence is preserved. For more complicated con¬ 
straints on which this procedure is not feasible, one 
can reject a proposed g' if it is not in the space of 
constrained networks. 

mi) In order to achieve a faster computation of s' = 
s(g') in 1(a) it is useful to store which motifs 
any given edge belongs to. Then, instead of re¬ 
computing all motifs of g' to compute s' = s(g '), 
one can calculate s' from the number of motifs de¬ 
stroyed and created when passing from g to g' using 
the procedure described in the note [u)]. 

iv ) A canonic ensemble with fixed /? is obtained by 
running loop 1. with a{s\ s) = ~P{s' — s). 

For a more formal description of the Wang-Landau al¬ 
gorithm and for further details we refer to Refs. [I, 2]. In 
Ref. [3j we provide an implementation of the algorithm 
described above for the case of triangles in undirected 
networks. 
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[2] D. Landau and K. Binder, A guide to Monte Carlo simu¬ 
lations in statistical physics (Cambridge University Press, 
2013). 

[3] A C++ code demonstrating the usage of Monte 
Carlo flat-histogram to sample constrained networks, 
https://dx.doi.org/10.5281/zenodo.30626. 
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